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Introduction 

Despite the fact that supercritical fluids occur both in 
nature and in industrial situations, the fundamentals of 
their behavior is poorly understood because supercrit- 
ical fluids combine the characteristics of both liquids 
and gases, and therefore their behavior is not intuitive. 
There are several specific reasons for the lack of under- 
standing: First, data from (mostly optical) measure- 
ments can be very misleading because regions of high 
density thus observed are frequently identified with liq- 
uids. A common misconception is that if in an ex- 
periment one can optically identify ‘drops’ and ‘liga- 
ments’, the observed fluid must be in a liquid state. 
This inference is incorrect because in fact optical mea- 
surements detect any large change (i.e. gradients) in 
density. Thus, the density ratio may be well below 
O(10 3 ) that characterizes its liquid/gas value, but the 
measurement will still identify a change in the index of 
refraction providing that the change is sudden (steep 
gradients). As shown by simulations of supercritical 
fluids [1], under certain conditions the density gradi- 
ents may remain large during the supercritical binary 
fluids mixing, thus making them optically identifiable. 
Therefore, there is no inconsistency between the opti- 
cal observation of high density regions and the fluids 
being in a supercritical state. A second misconception 
is that because a fluid has a liquid-like density, it is 
appropriate to model it as a liquid. However, such flu- 
ids may have liquid-like densities while their transport 
properties differ from those of a liquid. 

Considering that the critical pressure of most fuel 
hydrocarbons used in Diesel and gas turbine engines 
is in the range of 1.5 - 3 MPa, and the fact that the 
maximum pressure attained in these engines is about 6 
MPa, it is clear that the fuel in the combustion chamber 
will experience both subcritical and supercritical con- 
ditions. Studies of drop behavior over a wide range of 
pressures were performed in the past (Yang et al. [3], 
Delplanque and Sirignano [4|, Haldenwang et al. [5], 
and the review of Givler and Abraham [6]), however 
none of these studies identified the crucial differences 
between the subcritical and supercritical behavior. In 
tact, in two of these studies [3j. [5], it was found that 
the subcritical and supercritical behavior is similar as 
the drop diameter decreased according to the classical 


dMaw [7] over a wide range of pressures and drop di- 
ameters. 

The present study is devoted to the exploration of 
differences in fluid-behavior characteristics under sub- 
critical and supercritical conditions in the particular 
case of heptane fluid drops in nitrogen; these substances 
were selected because of the availability of experimental 
observations for model validation. 

Model equations 

The configuration studied is that of a single spherical 
drop in a medium with specified far field conditions. 
These far field values are identified by the subscript 
V and the location of the far field boundary, /?*(£), is 
calculated in a Lagrangian way to be that of null mass 
flux. 

The conservation equations are based upon Keiz- 
er’s fluctuation theory [8] which has the distinct ad- 
vantage of accounting for non-equilibrium processes. 
This formalism therefore leads to the most general fluid 
equations where the partial molar fluxes, 7*, and the 
heat flux, q , are related to thermodynamic quantities 
as follows: 

7 . = LtjVWnj), r = L„V0-^2 L qj V(0 h 

j J 

( 1 ) 

where 0 = 1 /(i^T), T is the temperature, is the 
universal gas constant, and fij are the chemical poten- 
tials. Here L x j are the Fick’s diffusion elements, L qq is 
the Fourier thermal diffusion element, L iq are the Soret 
diffusion, L qj are the Dufour diffusion elements, and the 
Onsager relations state that Lij — L JX and L xq — L 
Additionally, conservation of fluxes and mass in the sys- 
tem imply that m x J i — 0 and = 0 for 

j 6 [1, N] and j = q where m, are the molar masses and 
iV is the total number of species. 

Using the thermodynamic relationship 

iV-l 

' = J( Vjdp - hjdln T) + (£ a Djt dX t )/Xj (2) 

L 

where Uj is the partial volume, p is the pressure, X t is 


tin 1 molar traction, hj is the molar enthalpy, ami the 
mass diffusion factors are 

<*Dv ^liX&tJdXj =OXJdX J +X l dln yJdXj (3) 

one can calculate ./ t and If from 1 and 2. The con- 
servations equations were derived in detail in Harstad 
and Bellan [1] and results were obtained for supercrit- 
ical conditions. The emphasis is here on extending the 
calculations with the same model to the range of sub- 
critical conditions. 

Boundary conditions 

The detailed boundary conditions at r = R d have been 
derived in [1]; here we describe only the new aspects 
that enabled the calculations to be extended to sub- 
critical conditions. The jump conditions at the drop 
boundary are: mass balance; relationship between Rd 
and the emission flux F ems \ heat balance; balance of 
species 1 flux; and the nonequilibrium evaporation law. 
Additional equations at r = Rd are the momentum, 
and the equation of state for the mixture which is 
used twice (once on each side of the boundary). Thus 
there is a total of eight equations and nine unknowns: 
u b <^tb’^?b<Pb . Pb > Rd< Tb.Pb and F em , where the sub- 
script 6 denotes the drop boundary, p is the density, and 
superscripts L and G denote the inside and outside re- 
gions of the drop; variable ufc is obtained by integrating 
the drop continuity equation starting at r — 0. A ninth 
independent relationship exists only under subcritical 
conditions and is related to the existence of a surface, 
as discussed below. 

The indeterminacy of the boundary conditions for 
a fluid drop under supercritical conditions has already 
been discussed by Harstad and Bellan [1]. This is physi- 
cally understandable since there is no true surface, and 
thus there is an arbitrariness as to the choice of the 
boundary to follow. At least three choices are reason- 
able: One may follow the pure fluid boundary as was 
done by Harstad and Bellan [1]. Another possibility is 
to follow the initial boundary separating the two fluids, 
this being the choice in the present calculation. The 
third possibility is to follow the point of maximum den- 
sity gradient; although this is not the present choice, 
the point of maximum density gradient is calculated 
here & posteriori to indicate the location of the optically 
identified fluid drop. In contrast, under subcritical con- 
ditions the boundary to follow is the drop surface and 
the problem is fully determined. 

There are other important consequences of the ex- 
istence or lack of a surface at r = R^. For example, un- 
der strong evaporative conditions a mass fraction 'film’ 
layer exists inside the drop (9) and the thickness of this 
layer, 8y <C Ar“ where Ar~ is the distance from the 
surface to the first grid point inside the drop. A de- 
tailed analysis [2] shows that an effective mass diffusiv- 
ity D e /f can be defined with the consequence that a 


film layer exists when F t . ms > /;D,.///A/'~ . The valm 
of D e f f and that of an equivalent thermal conductivity 
Xp.fl were calculated under the quasi-steady assump- 
tion in [2] by finding two linear combinations of T and 
the mass fraction Y\ for which the transport matrix can 
be approximately diagonalized. In diagonal form, the 
characteristic length scales for diffusional transport of 
these two new variables are apparent, and this allows 
the definition of D t ff and A*// [2]. Previous calcu- 
lations [2] show that A e // > A and that D e ff < D. 
These definitions also allow the calculation of an ef- 
fective Lewis number Le e /f = A ef f /(nC p D e f f) once 
the values of the dependent variables are known. The 
quasi-steady assumption does not remove the generality 
of the estimate since the essence of the estimate is that 
of a characteristic length. One of the most important 
consequences of the mass fraction film layer existence 
is the direct relationship that exists between Y (Rd — e) 
and Y(Rd +£), where e max(<5y , 6x)\ it is this rela- 

tionship which provides the needed additional equation 
to fully determine the solution at the drop surface. This 
relationship can be formulated by considering the dif- 
ference AY X L = Y x (Rd - e) - Y x L (Rd - Ar~) where 
Y x (Rd — Ar“) represents the computational grid cen- 
ter value at the first adjacent position to the film layer 
inside the drop such that A r“ > 6y • Similarly one may 
define AT 1 = T L (Rd-e)-T L (Rd-&r-). The variable 
= exp (ipf - <pf), where <p is the fugacity, quantifies 
the Yj jump across the drop surface and can be cal- 
culated from the state equation. For example, under 
strict equilibrium evaporation (i.e. Fems = 0) condi- 
tions, = 1. For finite F erns and for a binary mixture 
system, its ratio to a reference state F rc /(£ 1? f 2 ) can be 
defined by 

Ferns = ^pFref (4) 

where 

Fref = SlOfav ~ *l) + B*(i 2 ~ £«J (5) 

Bj = a aj mjU Tj n, j = 1,2 ( 6 ) 

~ Ar~) - (7) 

A r-)]/[fl 2 y i i (/^ - A r-) + B x Y 2 L (R d - A r~)) 

where a aj is an accomodation coefficient, urj is the 
mean molecular velocity crossing the plane in one direc- 
tion, and n is the number of moles per unit volume [1]; 
consistently, F re /(1, 1) =0. A detailed analysis of the 
film layer yields then Xf(Rd-e),Xf(Rd + e).X$:(Rd- 
e) and X^iRdFe) in terms of e? and the £’s which pro- 
vides the additional relationship that allows closure of 
the system of equations at the drop boundary 

Since under supercritical conditions the concept of 
latent heat, and therefore of evaporation, is not ap- 
plicable, the above analysis does not hold. However, 
the film layer computational approach is still necessary 
if Pe g rtd > 0(1) in order to insure that all scales are 



n'solvod. Th<T<»for<\ the formalism of the tilm layer is 
retained for computational purposes even under super- 
critical conditions, although the layer no longer exists 
physically. Essentially, the solution in the supercritical 
regime has a diffusive character, whereas in the suberit- 
ical regime it has a diffusive-convective character where 
the convective part is introduced by the film layer and 
the evaporation. 

Results 

The present simulations are performed for an n-heptane 
drop in nitrogen because it is the set of binary sub- 
stances which is best documented experimentally. The 
equations of state have been calculated according to 
the procedure described in Harstad et al. [10], and the 
calculation of properties has been described in Harstad 
and Bellan [1]. The purpose of these simulations is to 
validate the model. 

The only data that can be used for comparisons 
is that obtained under evaporative rather than burning 
conditions, since in the last case the flame temperature 
that acts as the far field boundary is unknown. Fur- 
thermore, as shown below, it is only microgravity data 
that can be considered valid for these comparisons be- 
cause normal gravity data has unavoidable convective 
effects that are not modeled here. Additionally, since 
all high pressure microgravity drop evaporation experi- 
ments were performed with suspended drops, even these 
data are clearly not totally equivalent to our simulation 
results which are obtained for a free floating drop. 

To our knowledge, microgravity obtained data with 
C 7 H 16 drops evaporating in N 2 were reported only by 
Sato [11] and Nomura et al. [12]. In their experiments 
the 0.7-1 mm drops were suspended from a fiber of at 
least lOO^x diameter. The C 7 H 16 drop evaporation ex- 
periments of Chauveau et al. [13] were conducted only 
in normal gravity, whereas their reported micrograv- 
ity experiments were of burning drops. Therefore, our 
comparison focuses on the data of [11] and [12], while 
also considering for reference (see Table 1) the more re- 
cent normal gravity data of Morin et al. [14] for 1-1.5 
mm drops, instead of that of Chauveau et al. [13]. 

The simulations were performed for nominal initial 
conditions matching the experimental data: R% = 0.35 
mm except for the comparison with Sato’s [11] data 
which was performed for R° d = 0.5 mm, and = 300 
K. The far field conditions are located at R® = 4 mm 
where T e and p e are specified consistent with those of 
the experiments and V^ 0 = 0. The fluid drop is initially 
composed of pure heptane ( T c - 540.3 K,p c = 2.76 
MPa), while the surrounding is nitrogen ( T c = 126.2 
K. p c = 3.39 MPa); in order to avoid an initial unphys- 
ical discontinuity, a minute amount of heptane exists 
initially in the drop surroundings, its distribution van- 
ishing with increasing r. For the same reason, although 
the fluid drop temperature and drop-surroundings fluid 
composition are assumed initially uniform, a set of com- 


putational initial conditions (i.e. spatial profiles of tin* 
variables) are calculated for each simulation by satisfy- 
ing the nominal initial conditions at the domain bound- 
aries and the jump conditions at R,i . 

In all of the discussions below, subcritical’ and 
‘supercritical’ qualifications will be used with respect 
to the heptane critical point, and not with respect to 
the critical point of the mixture which varies according 
to the local composition. 

Determination of thermal diffusion factors tom high 
temperature data 

To proceed with the calculation, one must specify val- 
ues for the thermal diffusion factors, a, which can be 
defined either from the Irwing-Kirkwood (IK) or the 
Bearman- Kirkwood (BK) form of the heat flux [15] and 
ocbk ~ a iK ~ However, values of olqk are not 
well known for most substances, except at atmospheric 
conditions where they can be calculated from kinetic 
theory. Since we are here interested in calculations at 
considerably larger pressures, the question arises as how 
to calculate olqK' For this purpose, our premise is that 
if it can be shown that olik is very small, in fact it can 
be considered negligible with respect to ocbk ~ &ik, 
and then ocbk — which is only a function of ther- 
modynamic quantities [15]. Since is calculated from 
thermodynamics, this would provide an approximate 
value for ocbk f° r all (p, T) conditions where oci^jot* is 
very small having defined a* = max( p> x,x») | <*h |. The 
purpose of these high temperature data comparisons is 
to explore whether our premise that oljk/oc* is very 
small is correct. For heptane/nitrogen oc* = 0(10). 

Shown on Fig. 1 are (d/d 0 ) 2 plots from our simu- 
lations portraying Nomura et al.’s [12] experiments at 
high temperature (745 K) in the pressure range of 0.1 
- 2 MPa. In the numerical simulations, the location of 
the drop boundary is defined to be that of the maximum 
density gradient, for consistency with optical measure- 
ments. In agreement with well known theory [7], at 0.1 
MPa the liquid/gas interface is found to be precisely 
that of maximum density gradient. With increasing p 
the two locations still coincide for all simulations in the 
range 0.1-5 MPa investigated in this work, but the 
density gradient, although still substantial, decreases 
across the boundary as p increases. 

All but one of our simulations were conducted with 
ocjk = 0.01; the remaining simulation was conducted 
with ocqk = 0.01. Our results capture the 0.1 MPa data 
very well but display a somewhat earlier cpA aw behav- 
ior; it is unclear whether the non-coinciding part of the 
data and simulations fall within the experimental error 
since this error is not provided with the data. The data 
at 0.5 MPa is compared with results from simulations 
with both c*bk and a/# specified as 0.01. It is clear that 
the ocbk = 0.01 results fall short of agreement with the 
data, and in fact show a typical large increase in the 
evaporation time that was obtained with ocbk — 0.01 


ut i jf. lit*r pn\ssur<‘s .us well. In contrast. the e t / / x - 0.0 1 
results capture nonlinear portion of the curve very 
well with a small discrepancy in the total evaporation 
time. Simulations and data at 2 MPa agree only during 
the initial time, after which the simulations display the 
expected smooth variation consistent with drop heat- 
ing, whereas the data exhibit two discontinuities that 
can be explained only by the presence of the suspending 
fiber. Calculated slopes of the linear part of the curves, 
called the evaporation constant [7], K , are presented 
for comparison in Table 1 for the 0.1 and 0.5 MPa re- 
sults obtained with olik = 0.01. Despite the presence 
of the suspending fiber in the experiments, an excellent 
agreement exists between simulations and data. A sim- 
ilar comparison cannot be performed at 2 MPa since 
there is no evidence of linear behavior in the data. 

Confirmation of thermal diffusion /actors from interme- 
diate temperatures data 

Displayed in Fig. 2 are p = 2 MPa comparisons of sim- 
ulation results at 655 K for various values of a//c, one 
simulation where as k instead of a/*- is prescribed, and 
Nomura et al.’s [12] data at 656 K. The numerical pre- 
dictions are a very weak function of a//c in the range 
-0.6 - 0.6 and agree remarkably well with the data dur- 
ing the initial heat up period of the drop. Eventually, 
the data shows a faster evaporation than our simula- 
tions, although the lack of error bars in the data make 
it impossible to evaluate the extent of the disagreement. 
It is also difficult to evaluate the influence of the fiber 
(during the experiment) on the evaporation process. 
However, results with olbk = 0.01 clearly overestimate 
both the growth of the drop during the initial heat up 
time and the drop evaporation time. These results are 
consistent with those of Fig. 1. 

Additional comparisons between numerical predic- 
tions and data is portrayed in Fig. 3 where comparisons 
are made in the range 0.1-2 MPa between simulations 
at 655 K with a/# — 0.01, and data in the range 648 - 
669 K. The initial heating time is again very well repro- 
duced by the simulations, except that the predictions at 
0.1 MPa display again an earlier dMaw behavior. The 
evaporation time is very well reproduced at 0.1 MPa, 
and less well as the pressure increases. Since it is diffi- 
cult to quantify the influence of the suspending fiber as 
the pressure increases, we can qualify this comparison 
as very encouraging. 

Table 1 includes comparisons of K for this interme- 
diary temperature regime, and shows good to excellent 
agreement between data and predictions. 

This study indicates (see also below) that the value 
of ol[k /&* is indeed rather small and that qbk — a h, is 
correct. The assumption made in all calculations pre- 
sented below is that the value of c*/k is the same small 
value determined at high temperatures regardless of the 
(p<T) conditions, and thus that ciqk ~ a/j.This as- 
sumption might not be entirely valid, as in general a//c 


is a function of both /; and T This assumption and the 
fact that the data is from suspended drop experiments 
whereas our calculations are for free drops, might ax- 
plain the 15-20% discrepancies (see below and Table 1) 
between data and results from simulations. 

Comparison with data at low temperatures. 

The low temperature data of Nomura et al. [12] and Sato 
[11] (Sato’s data was approximated from his figure) is 
shown in Fig. 4 along with our numerical predictions 
at 445 K, 470 K and 495 K using olik ~ 0.01. The 
temperature range for Nomura et al.’s [12] data is 466 
- 493 K whereas Sato’s [11] data was obtained at 445 
K; the data in [12] is in the 0.1-5 MPa range, whereas 
that of [11] is at 2 MPa. The comparisons are very 
good at low p and deteriorate as p increases. The pre- 
dictions and data [12] agree remarkably well at 0.1 and 
0.5 MPa, whereas at 1 MPa the evaporation time is 
slightly overpredicted by the simulations. Nevertheless, 
the calculated and measured evaporation constant (Ta- 
ble 1) show very good agreement at all three pressures. 
The 2 MPa numerical results approximate the d 2 ex- 
perimental variation [11] fairly well, and the agreement 
in the value of K (Table 1) is excellent. At p = 5 
MPa, our simulation of a free drop shows an increased 
heating time, whereas the suspended drop in the exper- 
iment shows a decreased heating time with respect to 
the 0.1 MPa case. The difference between the experi- 
mental conditions and those of the simulations explains 
the disagreement in the heat up time, although the rate 
of regression of the largest gradient location is surpris- 
ingly well predicted. Since at 5 MPa the conditions are 
supercritical, there is no evaporation and the concept 
of evaporation constant is irrelevant, although compar- 
isons between the rates of regression are still meaning- 
ful. 

Conclusions 

A model of fluid behavior under both sub- and super- 
critical thermodynamic conditions has been discussed 
with particular emphasis on the different physics ac- 
cording to the initial conditions with respect to the crit- 
ical point. The model has been exercised for a fluid drop 
for which data are available for model validation. The 
drop is typically colder than its surroundings whose far 
field conditions are prescribed. In the subcritical regime 
and for large emission rates from the drop, there exists 
a film layer in the inner part of the drop surface and the 
solution of the equations has a convective-diffusive char- 
acter. In the supercritical regime, there is no material 
surface to follow, and this introduces an indeterminacy 
in the boundary conditions. To resolve this indetermi- 
nacy one must follow an arbitrary boundary of interest 
which is here that of the initial fluid drop. The solu- 
tion has then a pure diffusive character, and from this 
solution we calculate the location of the highest density 



gradient whit li wr identify with tin* optically observable 
fluid drop. 

The model was exorcised for a heptane drop in ni- 
trogen because of the existing data available for com- 
parison. Simulations obtained with this model were 
validated with microgravity experimental data for large 
drops over a wide range of temperatures and pressures. 
The large temperature data were used to determine the 
value of the thermal diffusion factor and further valida- 
tions were conducted with this fixed value. The agree- 
ment between predictions and d^data is excellent at at- 
mospheric pressure and becomes fair at supercritical 
pressure, whereas the rate of regression of the point of 
maximum density gradient is remarkably well predicted 
at ail pressures. The numerical predictions show that 
the traditional d 2 - law is obeyed only in the subcritical 
regime. As the pressure is increased, d 2 becomes non- 
monotonic with time, with a slope whose magnitude 
increases as a function of time. Thus, we initially iden- 
tify a heating period during which the drop size may 
increase, followed by a period during which the size is 
continuously reduced. The duration of the heat-up pe- 
riod increases with far field pressure. 
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Table 1: Maximum regression rate of the maximum density gradient location, K in mm 2 /s, obtained from the 
current model (ap), Nomura et al.’s microgravity experimental data (Nom), Sato’s microgravity and normal gravity 
experimental data (Sat), and Morin et ah’s normal-gravity data (Mor). The Nomura et al.’s and Morin et al.’s data 
were provided by the authors, and Sato’s values were read on their graph following the directions given in their 
paper. In the simulations Tj = 300K and d° = 0.7mm, while Nomura et al.’s dP was 0.6 - 0.8mm, Sato’s was 1mm, 


and Morin et al.’s was 1 - 1.5 mm. 



Fig. 1 High temperature comparisons. R% — 0.35 
mm; = 4 mm, y c ° = 0 and = 300 K. In the far 
field T e and p c are specified as in the experiments. Sim- 
ulations at T e = 745 K and p e : O.IMPa ; 0.5MPa, 

aix = 0.01 ; 0.5MPa, cxbk — 0 01 — ° — * 2MPa 

. Data: 741 K and O.IMPa ■; 749 K and O.SMPa 

A; 746K and 2MPa * 



Fig. 3 Intermediary temperature comparisons. 
= 0.35 mm; R® = 4 mm, Y e ° = 0 and b ~ 300 

K. Simulations at 655 K: O.IMPa ; 0.5MPa ; 

l MPa - ■ 2MPa . Data: 648 K and O.IMPa ■; 

655 K and 0.5MPa A; 669 K and IMPaT; 656 K and 
2MPa •. 



Fig. 2 Intermediary temperature comparisons at 
2MPa. i?3 = 0.35 mm; = 4 mm, Y° = 0 and 
T$ b = 300 K. Simulations at 655 K; aiK — 0.01 ; 

0.3 ; - 0.3 ; - 0.6 ;0.6 ; as * = 0.01 

— o — . Data at 656 K: ■. 



Fig. 4 Low temperature comparisons. R% - 0.35 
mm except at 445K where = 0.5 mm; R% = 4 mm, 
Y e ° = 0 and 7^ 6 = 300 K. Simulations at 470 K: 0.1 

MPa ; 0.5 MPa - - - ; 1 MPa ; at 445 K and 

2 MPa ; at 495 K and 5 MPa . Data: 471 

K and 0.1 MPa ■; 468 K and 0.5 MPa A; 466 K and 1 
MPa ▼; 445 K and 2 MPa 0; 452 K and 2MPa ►; 493 
K and 5 MPa •. 




